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Abstract. We study the long Josephson junction (LJJ) model which 
takes into account the second harmonic of the Fourier expansion of 
Josephson current. The sign of second harmonic is important for many 
p^ ' physical applications. The influence of the sign and value of the second 

harmonic on the magnetic flux distributions is investigated. At each step 
of numerical continuation in parameters of the model, the correspond- 
C I ing nonlinear boundary problem is solved on the basis of the continuous 

analog of Newton's method with the 4th order Numerov discretization 
scheme. New solutions which do not exist in the traditional model have 
Jh ^ been found. The influence of the second harmonic on stability of mag- 

y^i netic flux distributions for main solutions is investigated. 
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Ch I 1 Formulation of the problem 

Physical properties of magnetie flux in Josephson junctions ( JJs) play important 
role in the modern superconducting electronics. Tunnel SIS JJs arc known to be 
having the sinusoidal current phase relation. However, the decrease of the barrier 
transparency in the SIS JJs leads the deviations of the current-phase relation 
Q>^ I from the sinusoidal form [T]. 

VO ' We study the static magnetic flux distributions in the long JJs taking into 

lO I account the second harmonic in the Fourier-decomposition of the Josephson 

ly^ . current. This model is described by the double sine-Gordon equation (2SG) for 

^^ I magnetic flux distribution in the static regime [1],[3],[1]: 



o 



X 



— (^" + ai sin(^ + 02 sin2iy9 — 7 = 0, a;G(— /;/), (1) 

with the boundary conditions in the the following form 



H ■ ^'(±0 = /le. (2) 

Here and below the prime means a derivative with respect to the coordinate x. 
The magnitude 7 is the external current, I is the semilength of the junction, ai 
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and 02 are parameters corresponding the contribution of 1st and 2nd harmonic, 
respectively, h^ is external magnetic field. All the magnitudes are dimensionless. 

The sign of the second harmonic depends on a physical application under 
study. It is important, in particular, in junctions like SNINS and SFIFS, where N 
is a normal metal and F is a weak metallic ferromagnet [5] . Interesting properties 
of long Josephson junctions with an arbitrarily strong amplitude of the second 
harmonic in current phase relation were considered in [6]. We investigate the 
existence and stability magnetic flux distributions in dependence on the second 
harmonic contribution in both cases of negative and positive sign. 

Stability analysis of (f{x,p) is based on numerical solution of the correspond- 
ing Sturm-Liouville problem |7l8j : 

-^" + q{x)ij = \lP, V'(±0 = (3) 

with a potential q{x) — ai cos ^p + 2a2 cos 2ip. 

The minimal eigenvalue Ao(p) > corresponds the stable solution. In case 
\{p) < solution ip{x,p) is unstable. The case Xa{p) = indicates the bifurca- 
tion with respect to one of parameters p = (?, ai, 02, ft-e, 7)- 



2 Numerical scheme 

For numerical solution of the boundary problem ([II),© we apply an iteration 
algorithm based on the continuous analog of Newton's method (CANM) [5]. 
Let an initial approximation for (po{x) be given. At fc*'' step (fc = 1, 2, , . . .) we 
calculate: 

1. Iteration correction Wk{x) by solving linearized boundary problem 

-Wk+ qk-i{x)wk = <^fc_i - fk-i{x) , (4) 

w'k{-l) = -ip'k_,{-l) + K,, (5) 

w'kil) ^-f'k-lil) + he, (6) 

where f{x) = ai sin ip + a2 sin 2ip — 7. 

2. Next approximation 

ipk{x) == ipk-lix) +TkWk 

where parameter r^ is calculated by the Ermakov-Kalitkin formula [9] . 

Further in order to simplify notations the iteration indices are omitted. 

We introduce the grid M^ = {xi = —l + {i — l)h,i = 1,N,xn = l,h = 
21/{N—1)}. Numerov's discrete approximation [TO] of (H])-® yields the following 
linear algebraic system with the three-diagonal structure at i = 3, A^ — 2: 

—25wi + 48w2 — 36w3 + I6W4 — iw^ = 12/i(/ie — fi) 

a2Wi + 62W2 + C2W3 + d2W4 + e2W5 = r2 
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16.2809 



Table 1. Values of function ip and quantities tr ([7]) in some points of the interval 
[-/; /] for solution of kind <?^ at 21 = 10, 7 = 0, he = 0, ai = 1, 02 = 0. 



aiWi^i + hiWi + CiWi+i = ri, i = 3,N ~2, 

Gn-IWn + bN-lWN-l + Cn-IWn-2 + dN-lWN-3 + ^N-lWN-i = ''TV-l 
2bWN — 48u'Ar_i + 36'U;Ar_2 ~ 16wAr-3 + 3u)jv_4 ~ 12h{he, — ^'n) 

where the coefficients are determined by the following way 

7/j2 5^2 ^2 ^2 

02 = 1 , 62 = -2 ^ 92 , C2 = 1 + — 93 , d-i^ — - 94 , 62 = — 95 , 

r2 = ^ (I4/2 - 5/3 + 4/4 - /s) - ^ (14^ 'i - 5^ 'i + ^^l~^ 'D , 



12 



12 



/l2 



he 



flz = 1 - Y^ qt-i , h.i = -2 ^ g, , Ci = 1 - Y^ gi+i , i = 3, TV - 2 , 



^. = Y^ ih+i + 10/. + ./n) - ^ iVi+i + 10^.r + V'TJ , * = 3, TV - 2 . 

aw-1 — 1 , ow-i — —2 — gw-i , cn-1 — 1 + —77 'Zjv-2 , 

O 12 






h^ h^ 

-Z- qN-3 , eAr_i = — gAr_4 , 



rAr_i = — {-JN-i + 4/Ar_3 - 5/Ar_2 + 14/jV-l)- 



-- (-^" + 4^^_3 - 5(^'^_2 + 14^^-1) , 

where ipi = (p{xi), q, = q{x.j), f, = f{xi). 

In order to test the accuracy order of the above numerical scheme we perform 
the calculations of ([T]) , ([2]) at the sequence of uniform grids with steps h, h/2 and 
h/A {h ~ 0.15625). The results for solutions of the kind <?^ arc presented in the 
table [TJ It is seen, the quantities a calculated by formula 



, X ^h{x,) - iph/2iX;) . 

a(Xi} = ^— , 1 = 1,2,. 

fh/2[X^) -(Ph/4[X,) 



.,iV, 



(7) 
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are close to 2^ that corresponds the theoretical accuracy order 0{h^) of Nu- 
merov's approximation. 

The StLP ([3]) is approximated by the three-point finite-difference second 
order formulas |11| . The resulting algebraic eigenvalue problem is solved numer- 
ically with help of a standard subroutine [12] • 

3 Numerical results and conclusions 



Trivial solutions. In the "traditional" case 02 = two trivial solutions (p = 
and if =: TT (denoted by A/o and M^^ respectively) are known at 7 = and he = 0. 
Accounting of the second harmonic 02 sin 2ip leads to appearing two additional 
solutions (f = ±arccos(— 01/202) (denoted as M±ac)- The corresponding Aq as 
functions of 2SG-equation coefficients have the form Aq [Mq] = oi +2a2, Aq [M-^] = 
—ai + 2a2 and Ao[Af±oc] = {of — 4a|)/2a2. The exponential stability of these 
constant solutions (CS) is determined by the signs of the parameters ai and 02 
and by its ratio 01/02. 

The full energy associated with the distribution of ip{x) is calculated by the 
formula [5] 



F{p) 



1 , 

2^ 



+ 1 - qix) - j(p 



dx — heAip. 



The full energy behavior in dependence on 02 for considered distributions in 
the junction at /le = 0, 7 = 0, ai ~ 1, 21 ^ 10 is shown in Fig. [TJ 

Stability properties of trivial solutions have been investigated in |13j . 
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Fig. 1. Full energy in dependence on 02 £ [—1; 1] at /ig = 0, 7 = and 21 = 10 
forCS and^i. 



Fluxon solutions. The fluxons play an important role in the JJ physics. 
At small external fields h^ such distributions are fluxon <?^, antifluxon <P~^ and 
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l:ai = 
"^'^^ 2:^2 = 0.2 

3: (32 = 0.5 
4: fl, = 0.7 




Fig. 2. Dependence Ao(/ie) for <1>^ at increase of 02 £ [0; 1] and ai = 1, 2Z = 10, 

7 = 0. 

their bound states $^$~^ and (1>~^<P^ . As external magnetic field h^ is growing, 
more complicated stable fluxon and bound states appear: <?^" and ^='="^T« 
(n = 1,2,3,...). 

The energy of one-fluxon distribution <P^ limits to unit F{a2 — ?► 0) ^ 1 which 
corresponds to an energy of a single fluxon <P]^ in a traditional "infinite" junction 
model at oi = 1, 02 = 0. With change of 02 the number of fluxons ,8^ 



Nip) 



2/7r 



(p{x) dx, 



corresponding to the distribution <P^ is conserved i.e. dN/da2 = 0. Here we have 
N[<P^] = 1. 

Let us discuss the main features of the dependence Ao(/ie) for one-fluxon 
state (l>^ in two intervals: 02 G [0, 1] and 02 G [—1,0]. 

The change of the curve Xo{he) for one-fluxon state <P^ when the parameter 
02 increases in the interval 02 S [0; 1] is shown in Fig. [2] When he = 0, the 
state <?^ remains unstable. With increase in 02, Aq increases monotonically and 
tends to zero. With increase in magnetic field this solution becomes stable. The 
bifurcation point moves to the left with increase of parameter in the interval 
a2 G [0; 0.7]. At 02 > 0.7 it goes to the right again. It follows a stability interval 
ending at her ~ 2. The second bifurcation point also moves to the left at 02 G 
[0;0.5]. When 02 is increased in 02 G [0.5; 1], the bifurcation value her is also 
increased. 

In the interval 02 G [—1,0] we observe the following. When 02 increases in 
(—0.5; 0], the curve Xo{he) for #^ moves to the right (Fig.|3]). At 02 < —0.5 (case 
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Fig. 3. Xo{he) in dependence on 02 G [-1; 0] for (p^ and <?^*, at oi = 1, 2Z = 10, 

7 = 0. 



a2 = —0.7 in Fig. |3]), the curve corresponding to the stable solution <l>^ has two 
separate branches that are intersected at 02 ~ —0.8). Here we observe a region 
along he, where two different stable one-fluxon solutions (denoted by ^^ and 
#^*) coexist, see Fig. S) 

Thus, we considered both positive and negative contributions of the second 
harmonic in 2GS equation. It is shown that its accounting leads appearing new 
constant solutions and changes the stability properties of the fluxon solutions. 
Coexisting of two stable one-fluxon solutions requires further analysis and phys- 
ical interpretation. 
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